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Abstract 

The organization of the axonal cytoskeleton is a key determinant of the normal function of an 
axon, which is a long thin projection away from a neuron. Under normal conditions two axonal 
cytoskeletal polymers microtubules and neurofilaments align longitudinally in axons and are 
interspersed in axonal cross-sections. However, in many neurotoxic and neurodegenerative 
disorders, microtubules and neurofilaments segregate apart from each other, with microtubules 
and membranous organelles clustered centrally and neurofilaments displaced to the periphery. 
This striking segregation precedes abnormal and excessive neurofilament accumulation in these 
diseases, which in turn leads to focal axonal swellings. While neurofilament accumulation 
suggests the impairment of neurofilament transport along axons, the underlying mechanism of 
their segregation from microtubules remains poorly understood for over 30 years.To address this 
question, we developed a stochastic multiscale model for the cross-sectional distribution of 
microtubules and neurofilaments in axons. The model describes microtubules, neurofilaments 
and organelles as interacting particles in a 2D cross-section, and is built upon molecular 
processes that occur on a time scale of seconds or shorter. It incorporates the longitudinal 
transport of neurofilaments and organelles through this domain by allowing stochastic arrival 
and departure of these cargoes, and integrates the dynamic interactions of these cargoes with 
microtubules mediated by molecular motors. Simulations of the model demonstrate that 
organelles can pull nearby microtubules together, and in the absence of neurofilament transport, 
this mechanism gradually segregates microtubules from neurofilaments on a time scale of hours, 
similar to that observed in toxic neuropathies. This suggests that the microtubule-neurofilament 
segregation is simply a consequence of the selective impairment of neurofilament transport. The 
model generates the experimentally testable prediction that the rate and extent of segregation 
will be dependent on the sizes of the moving organelles as well as the density of their traffic. 


Author Summary 

The shape and function of axons is dependent on a dynamic system of microscopic intracellular 
protein polymers (microtubules, neurofilaments and microfilaments) that comprise the axonal 
cytoskeleton. Neurofilaments are cargoes of intracellular transport that move along microtubule 
tracks, and they accumulate abnormally in axons in many neurotoxic and neurodegenerative 
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disorders. Intriguingly, it has been reported that neurofilaments and microtubules, which are 
normally interspersed in axonal cross-sections, often segregate apart from each other in these 
disorders, which is something that is never observed in healthy axons. Here we describe a 
stochastic multiscale computational model that explains the mechanism of this striking 
segregation and offers insights into the mechanism of neurofilament accumulation in disease. 


Introduction 


Axons are long slender projections of nerve cells that permit fast and specific electrical 
communication with other cells over long distances. The ability of nerve cells to extend and 
maintain these processes is critically dependent on the cytoskeleton, which is a dynamic scaffold 
of microscopic protein polymers found in the cytoplasm of all eukaryotic cells. The axonal 
cytoskeleton comprises microtubules, intermediate filaments called neurofilaments, and 
microfilaments. Microtubules and neurofilaments are both long polymers that align in parallel 
along the long axis of the axon, forming a continuous overlapping array that extends from the 
cell body to the axon tip Qi- Microtubules are stiff hollow cylindrical structures about 25 nm 
in diameter that serve as tracks for the long-range bidirectional movement of intracellular 
membranous organelles and macromolecular cargo complexes. In axons, this movement is 
known as axonal transport Q. The cargoes of axonal transport are conveyed by microtubule 
motor proteins: kinesins in the anterograde direction (towards the axon tip), and dyneins in the 
retrograde direction (towards the cell body) Q. Neurofilaments, which are the intermediate 
filaments of nerve cells, are flexible rope-like polymers that measure about 10 nm in 
diameter Q. These polymers function as space-filling structures that expand axonal 
cross-sectional area, thereby maximizing the rate of propagation of the nerve impulse ||^[^. In 
large axons, neurofilaments are the single most abundant structure and occupy most of the 
axonal volume Q. Mutant animals that lack neurofilaments develop smaller caliber axons and 
exhibit delayed conduction velocities SOI- 

In addition to their structural role in axons, neurofilaments are also cargoes of axonal 


transport, moving along microtubule tracks powered by kinesin and dynein motors 113 ■ 16 1261. 
The filaments move at rates similar to membranous organelles but the movements are less 
frequent, resulting in a “stop and go” motile behavior characterized by short bouts of movement 
interrupted by prolonged pauses on a time scale of seconds or shorter IIZ10- The net result is 
an average rate of transport that is much slower than that for many other cargoes. 

Neurofilaments has been observed to accumulate abnormally in axons in many 
neurodegenerative diseases including amyotrophic lateral sclerosis, hereditary spastic paraplegia, 
giant axonal neuropathy and Charcot-Marie-Tooth disease (also known as hereditary distal 
motor and sensory neuropathy) ||^ 19 ■^, and also in many toxic neuropathies |24-^. In 
extreme cases, these accumulations can lead to giant balloon-like axonal swellings | |29f[34) . 
These accumulations are thought to be caused by alterations in neurofilament transport, but the 
mechanism is not understood | |T5| - |40l . 

In healthy axons, microtubules and neurofilaments align along the length of an axon and are 
interspersed in axonal cross-sections 11 41 42} 143JI, with microtubules often forming small 
clusters in the vicinity of membranous organelles ||8 44 451. However, in many toxic and 
neurodegenerative disorders these polymers segregate, with microtubules and membranous 
organelles typically clustered in the center of the axon, and neurofilaments displaced to the 
periphery (Fig. [^. This striking cytoskeletal reorganization, which is never observed in healthy 
axons, has been reported in neurodegenerative disorders as diverse as giant axonal 
neuropathy |46-^ and Charcot-Marie-Tooth disease |34 as well as in neurotoxic 
neuropathies induced by exposure to agents as diverse as 2,5-hexanedione and 
3,3’-iminodiproprionitrile (IDPN) |24 5T-|57} 1311, aluminum | [58) , carbon disulfide |59 60|, 
estramustine phosphate | [6T[ , 1,2-diacetylbenzene ^ 2 \ and 1,2,4-triethylbenzene | |63j , and in a 
transgenic mouse expressing a mutant neurofilament protein |64|. However, the mechanism of 
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this segregation and its relationship to the neurofilament accumulation that also occurs in these 
different conditions is not known. 



Figure 1. Cross-sectional distributions of microtubules and neurofilaments via electron 
microscopy. (A) Schematic drawing that illustrates the normal distribution of microtubules 
(large black dots), neurofilaments (small grey dots) and organelles (cyan disks) in untreated 
axons. (B) Schematic drawing that illustrates the segregated components in IDPN-treated axons. 
(C) In normal axons microtubules, neurofilaments and organelles are interspersed. Big dots: 
microtubules; small dots: neurofilaments; circular objects: organelles. (D) In IDPN-treated 
axons these components segregate in a core-tube pattern, with microtubules and organelles 
typically in the center and neurofilaments displaced to the periphery. The scale bars are 1 /im. 
Reproduced from Papasozomenos et al, 1981. Originally published in Journal of Cell Biology, 
91:866-871. 


Microtubule-neurofilament segregation has been studied most extensively for IDPN and 

2.5- hexanedione. IDPN is a compound closely related to the naturally occurring food poison 
3-aminopropionitrile which causes the neurological disorder lathyrism |65 ■^, and 

2.5- hexanedione is a metabolite of the industrial solvent hexane. The mechanism of toxicity is 
not known, but it is thought to involve chemical modification of neurofilaments, which 
presumably disrupts their normal interactions with microtubules in some way |25 28 69 ■ 74|. 
Systemic administration of IDPN or 2,5-hexanedione to rats causes selective impairment of 
neurofilament transport |75 -7^, focal accumulations of axonal neurofilaments leading to axon 
enlargement, and neurological defects similar to amyotrophic lateral sclerosis (ALS) in 
humans | [80|j^ . Injection of IDPN or 2,5-hexanedione into peripheral nerves results in local 
microtubule-neurofilament segregation within just a few hours, preceding the accumulation of 
neurofilaments by hours or days 151 52 56 84 131) . This segregation does not appear to affect 
the axonal transport of membranous organelles, which continue to interact with and move along 
these tracks in spite of their clustering. Moreover, in the case of IDPN the segregation has been 
shown to be reversible p4||131|, as has the impairment of neurofilament transport | |85) . In | 
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a single injection of IDPN into rat sciatic nerves resulted in segregation in axons at the injection 
site within a few hours, but the segregation disappeared in about a day. In 11311, a single 
injection of IDPN into the body cavity of rats resulted in segregation within the axons of the 
sciatic nerve after 4 days, and this disappeared after six weeks. Thus the 
microtubule-neurofilament segregation caused by IDPN and 2,5-hexanedione is fast, local and 
reversible. 

Though the segregation of microtubules and neurofilaments in axons was first described 
more than 30 years ago, the underlying mechanisms are still poorly understood. Given that 
neurofilaments move along microtubule tracks and that microtubule-neurofilament segregation 
precedes neurofilament accumulation and axonal enlargement in rodent models, it is attractive to 
speculate that the segregation reflects an uncoupling of neurofilaments from their transport 
machinery CD- However, the mechanism by which such an uncoupling at the molecular level 
might generate polymer segregation at the population level remains unclear. 

To address these questions, we have developed a stochastic multiscale model for the 
cross-sectional organization of microtubules and neurofilaments in axons. The model describes 
microtubules, neurofilaments, and organelles as interacting particles that move in a 2D domain 
representing a cross-section of an axon, and incorporates axonal transports of neurofilaments 
and organelles, as well as volume exclusion and Brownian motion of all the particles. 
Neurofilaments and organelles dynamically bind to and unbind from nearby microtubules 
through molecular motors, and the motor cross-bridges are modeled as elastic springs. The 
longitudinal movement of neurofilaments and organelles along axons is modeled by stochastic 
addition and removal of these cargoes. The multiscale nature of the model lies in that it is built 
upon molecular processes that occur on a time scale of seconds or fractions of a second, and 
addresses segregation phenomena of two populations of polymers that occur on a time scale of 
hours to days. 

Simulations of the model demonstrate that if we block neurofilament transport by preventing 
neurofilament from binding to microtubules, then organelles pull nearby microtubules together 
and gradually segregate them from neurofilaments, on the same time course as observed in toxic 
neuropathies; while if we restore neurofilament transport, then microtubules and neurofilaments 
start to remix until their spatial distribution returns to normal. The model further predicts that (1) 
during the segregation process, microtubules first form small clusters, small clusters merge into 
bigger clusters, and eventually a single cluster forms close to the center of the domain; (2) in the 
absence of neurofilament transport, larger organelles are more effective in causing complete 
cytoskeletal segregation than small organelles with the same density. 


Model 

The stochastic multiscale model 

In our model, microtubules, neurofilaments and organelles are described as individual particles 
that move in a circular domain D with fixed radius Rq, representing a cross-section of an axon. 
Microtubules and neurofilaments are rod-like polymers that align along the length of axons, thus 
they are treated as nondeformable disks in D (Fig. |^), with center positions denoted by 
yf) and radii by r^. Here k = M or N is the index for particle type: M for 
microtubule, N for neurofilament; and i with 0 < i < is the index for the A:-type particle 
where is the total number of fc-type particles. The radii of microtubules and neurofilaments 
are constant, with = 12.5 nm and = 5 nm. Organelles in axons have different sizes and 
shapes, and their cross-sectional geometry depends on their position relative to the cross-section 
(Fig. 1^). In this model, we took organelles as spindle-shaped objects and, for simplicity, we did 
not consider possible shape changes (Fig. |^). Therefore the organelles exist as non-deformable 
disks in D, and as an organelle crosses D, its cross-sectional radius, r^, varies according to its 
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position, z2, relative to D, 


= 6 1 - 




—a < zP < a 


( 1 ) 


Here a is half of the organelle length, 6 is its maximum cross-sectional radius, z^ is the distance 
of its center to D, and the index “O” stands for organelle. By varying the parameters a and 6, we 
can vary the overall dimension of the organelles (Fig. |^). 

We examined three key molecular mechanisms that contribute to the cross-sectional 
distribution of microtubules and neurofilaments: slow axonal transport of neurofilaments, fast 
axonal transport of organelles, and volume exclusion of all the particles. In the following 
sections we describe in detail how these mechanisms were incorporated into our model. We 
denote the unit vector pointing from to x* by and the surface distance between the i-th 
particle of fc-type and j-th particle of ^-type by d^, given as. 


d'lj = |xf -x'l -rf -r'. 


( 2 ) 


Mechanism 1: Slow axonal transport of neurofilaments. 


Neurofilaments interact with molecular motors (kinesin and dynein) which move these polymers 
along microtubules either anterogradely or retrogradely 1^^^. The movements are fast 
but infrequent because the filaments spend most of their time pausing, which results in a slow 
average rate of transport 188 1126| . This longitudinal movement of neurofilaments along 
microtubules can change the spatial distribution of these polymers in axonal cross-sections. For 
example, if a neurofilament moves into the cross-section along a microtubule in a small cluster 
of microtubules it can displace one or more of the adjacent microtubules, dispersing the cluster. 
Alternatively, if a neurofilament between two microtubules moves out of the cross-section, then 
the two microtubules are able to diffuse closer together. 




Figure 2. Model setup. (A) The model geometry. The computational domain D represents a 
cross-section of an axon. The black circle is the domain boundary representing the axon 
membrane. Small grey dots are neurofilaments (NF), large black dots are microtubules (MT) 
and cyan filled disks are organelles (Org). (B) The relation of D to the whole axon. Thin grey 
lines are neurofilaments, thick black lines are microtubules, the cyan body is an organelle, red 
triangles represent molecular motors that move microtubules and organelles along microtubule 
tracks. (C) The shape of organelles considered in this model. The cross-sectional radius of an 
organelle in D depends on its position relative to D. 

Based on the above considerations, we modeled neurofilament transport in the following 
way. A neurofilament in D can bind to a microtubule within a binding radius with rates 
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A neurofilament bound to a nearby microtubule can unbind with rate k^ff, or move away and 
exit D with rate We conserved neurofilament number in D by replacing each departing 
neurofilament with a new entering neurofilament, placed at a distance i?;, from a randomly 
selected microtubule. Neurofilaments and microtubules are long polymers that are aligned in 
parallel along the long axis of the axon so a new neurofilament can only enter D if a 
neurofilament was present at that location in the adjacent plane at the preceding time point. 
Thus, to prevent the entry of a new neurofilament in a region of D that is lacking other 
neurofilaments (specifically, this would be encountered when simulating the remixing of 
neurofilaments and microtubules after segregation), we only permitted the entry of a new 
neurofilament next to a microtubule that was already within a radius Rb of another 
neurofilament. We did not differentiate the direction of neurofilament movement along the axon 
because anterograde and retrograde movements have similar contributions to the distribution of 
neurofilaments in D. If a neurofilament is bound to a microtubule, they interact through the 
following elastic spring forces. 


r^MN 


r^NM 

^j,i 


_ N^MNMN 

Uj^j tJjji 


(3) 


Here and are the forces acting on the i-th microtubule by the j-th neurofilament 

and vice versa, and is the spring constant. If is bigger than the binding radius Rb, then 
there is no spring force between the neurofilament-microtubule pair. 

We assumed that each neurofilament could engage with only one microtubule at a time. The 
rationale for this is as follows. The on-rate and off-rate for neurofilament binding to 
microtubules is estimated to be 10“^/s and 6.5 x 10“^/s based on previous experiments. Thus 
a neurofilament within the binding radius of a microtubule would spend, on average, 

1/(1 -t- 6.5) « 0.13 of its time engaged with that microtubule, and the chance for one 
neurofilament within the binding radius of two microtubules to bind both simultaneously would 
be 0.13^ « 0.017. Since in reality no neurofilament would remain within the binding radius of 
two microtubules at all times, the actual probability is even lower. Thus the chance for one 
neurofilament to interact with multiple microtubules simultaneously is negligible and we neglect 
it in our model. 


Mechanism 2: Fast axonal transport of organelles. 


Like neurofilaments, membraneous organelles are also conveyed anterogradely or retrogradely 
along microtubule tracks by kinesin and dynein motors. However, these cargoes tend to spend 
much less time pausing, resulting in a much faster average rate of transport. Due to their large 
size, the movement of organelles can cause significant fluctuations of the microtubule and 
neurofilament organization by displacing these polymers laterally. These cargoes can bind 


multiple motors 189 901 and, due to their large size, they can readily interact with multiple 
microtubules even if those microtubules are not close to each other | [45l[M| . As an organelle 
moves along several microtubules, it can pull them closer together, similar to a “zipper”. This 
speculation is supported by in vivo data that demonstrate organelles being surrounded by 
multiple microtubules in close proximity ||8 44 45 911, and in vitro experiments | |9^ which 
show that motors bound to spherical cargoes can pull on multiple microtubules and align them. 

Based on the above considerations, we modeled organelle movement in D in the following 
way. Organelles enter D randomly with rate kf^, moving along randomly chosen microtubules, 
and move persistently until they leave D completely. Thus each organelle is present in D for a 
time period that equals its length (2a) divided by its speed and 


where t is the time the organelle has been present in D. As the organelle moves from one side of 
D to the other, its cross-sectional radius rf first increases from 0, reaches its maximum when it 
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is halfway through, and then decreases to 0, which is given by Eqn Q . While present in D, an 
organelle can bind stochastically to an available microtubule within a binding radius Ri, with 
rate and unbind with rate j. If an organelle and a microtubule are bound, they interact 
through the linear spring force. 


r<MF 


r^FM _ jMF^MF 


(5) 


Here and are the forces acting on the i-th microtubule by the j-th organelle and 

vice versa, and k® is the effective spring constant which represents the action of possibly 
multiple motors. 


Mechanism 3: Volume exclusion. 


In addition to the active movement of neurofilaments and organelles and their interactions with 
microtubules through molecular motors, all the particles in the system interact through forces of 
volume exclusion. 

Neurofilaments have sidearms which are highly-charged unstructured polypeptide domains. 
These sidearms project outward from the filament core to form an entropic brush that defines a 
zone of exclusion around the polymer via long-range repulsive forces |93 ■^, maximizing the 
space-filling properties of these cytoskeletal elements. Microtubule associated proteins such as 
tau also have highly charged long polypeptide domains that can have a similar volume-excluding 
effect 193 9^ 100| . Based on these biological considerations, we modeled volume exclusion of 
neurofilaments, microtubules and organelles through the following pairwise repulsions. 


T>kl _ 
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if 


if dfi > Lr 


(6) 


Here is the force acting on the i-th particle of fc-type by the j-th particle of Z-type, where 
k, I = M, N or O, 1 < i < and 1 < j < For example, is the force acting on the 
i-th microtubule by the j-th neurofilament. Here Lr is the maximum interaction distance; 
specifies the magnitude of the force; and the negative sign preceding indicates that the force 
is repulsive. We note that this force goes to infinity as the surfaces of two particles approach 
each other and remains zero if the distance between two particles is larger than Lr- The 
functional form of the force is similar to those used in 1 101[ 122| for neurofilament repulsions 
and matches recent experimental data | |103| . 

To keep all the particles inside the domain, we modeled volume exclusion of the particles 
with the domain boundary in a similar way. The force acting on the i-th particle of fc-type by the 
axonal membrane is given by 


r: 


kB 


- (hrI dfl^ - l) if dfl^ < Lr 

0 if > Lr. 


(7) 


Here the index B stands for “boundary”, d^^ = Rq — |x^| — and is the unit vector 
pointing from the center of the domain to x^. 

Microtubules, neurofilaments, and organelles can also interact with each other 
hydrodynamically through the axoplasm. Organelle movement can cause significant flow of the 
axoplasm near their surfaces and displace nearby microtubules and neurofilaments. As an 
organelle pushes into D, its radius increases and it pushes nearby fluid and particles away from 
itself; as it moves away from D, instead of leaving void behind it, it creates negative pressure 
which draws the axoplasm to flow back and fill the space. The hydrodynamic effect due to the 
movement of microtubules and neurofilaments is presumable smaller given their constant and 
smaller size in cross-section. In this model, we do not model the hydrodynamic interactions 
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among these particles explicitly, but include this effect by adjusting the force prefactors 
associated with organelles. Specifically, when an organelle push into the domain, we double 
and to take into account the contribution of the fluid flow it creates. 


Model equations. 

The movements of microtubules, neurofilaments, and organelles in axons are viscous-dominated 
and thus inertia can be neglected. Under this simplification, we have the following system of 
stochastic differential equations 


dxf = dt + cTk dWj^, 1 < i < n'^, k = M,N,F. 


(8) 


Here F^ is the sum of all applied forces on that particle specified in P]),©, ® and Q. For 
example, Ff = E, + E, Rff + + E, + Rf ^The constant 

denotes the drag coefficient of the A:-type particle. Finally Wf are independent 2D Wiener 
processes modeling the random motion of these particles, and the amplitude is given by 
CTfc = \/2Dk, where Dk is the diffusion coefficient of the particle calculated by the Einstein 
relation. 


Parameter estimation and simulation algorithm 

The parameters used in our model are physical, and thus they are all measurable. Most of them 
have already been measured or there exist experiments that can be used to estimate them. Table 
[^summarizes all the parameter values, and the detailed methods to obtain these parameters are 
given in the SI. The units of these parameters reflect the time scales for the molecular processes 
integrated into the model, which are seconds or fractions of a second. 

To solve the model numerically, we treated the binding and unbinding, arrival and departure 
of cargoes explicitly at discrete time steps, and integrated the model system (j^ using the explicit 
Euler’s method. Because ak, k = M, N, C are constant in time, the numerical integrator has 
strong order 1.0 fMT) . We chose a time step h much smaller than all the time scales involved in 
Mechanisms 1-3. For the simulations of segregation and remixing over hours to a day, we used 
h = 1/50 sec if there was no organelle in D, and h = 1/1600 sec otherwise in order to deal 
with the stiffness of the equations introduced by the pushing of organelles when they move into 
D. The detailed simulation algorithm is included in the SI. The computational tool is written in 
C-H-. 


Results 


The organization of neurofilaments in normal axons. 

Morphometric studies suggest that neurofilaments are spaced randomly in axonal cross-sections 
when packed at low densities, but as the density increases they start to experience the 
volume-exclusionary repulsive forces of their neighbors and assume a less random distribution 
characterized by a more even neurofilament spacing |41 1122 143| . In this section we 
demonstrate that the neurofilament distribution generated using our model agrees well with 
these experimental data. 

Different methods have been used to characterize neurofilament distribution in axonal 
cross-sections. Kumar et al | |122| | used the radial distribution function (RDF) (also known as the 
pairwise correlation function). The RDF, denoted as g{r) describes how density varies as a 
function of distance from a reference particle. For particles that move randomly and completely 
independently, g{r) is a constant value of 1; while for crystalline structures g{r) forms peaks at 
precisely defined intervals. For neurofilaments in axons the shape of g{r) typically lies between 
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Table 1. Model Parameter Values 


Parameter 

Description 

Values 

/w 

Radius of neurofilament backbone 

5 nm 

J.M 

Radius of microtubule backbone 

12.5 nm 

Rb 

capturing radius for microtubule-cargoes 
active binding 

80 nm 

lTV 

^on 

rate for neurofilament binding 

1.0 X 10-Vs 

uN 

'^off 

rate for neurofilament unbinding 

6.5 X 10-Vs 

uN 

^out 

rate for neurofilament departure 

0.1 s-i 

^on 

rate for organelle binding 

2/s 

^in 

rate for organelle unbinding 

2/s 

rate for organelle passage 

0.105/s 

sO 

speed of organelle movement along mi¬ 
crotubules 

1 /rm/s 

Ij j' 

characteristic repulsion distance 

121.2 nm 


repulsion scale ( = 

0.5 pN 


effective spring constant for microtubule- 
neurofilament binding 

0.18pN/nm 

kO 

effective spring constant for organelle- 
microtubule binding 

0.9pN/nm 


drag coefficient of neurofilaments 

73.5 pN • s/^m 


drag coefficient of microtubules 

512 pN • sVm 


drag coefficient of organelles 

40.3pN • s//rm 

Dn 

diffusion coefficient of neurofilaments 

5.59 X 10-VwVs 

Dm 

diffusion coefficient of microtubules 

8.02 X lO-S/rmVs 

Dp 

diffusion coefficient of organelles 

1.02 X lO—^/rm^/s 


E.E.: estimated from experiments; see SI for detailed information. 


Notes and Refs 
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these two extremes, increasing sharply from 0 and forming a peak around 30 — 50 nm | 
Another method used often is to calculate the occupancy probability distribution (OPD), which 
is the distribution for the number of particles within an observation window of a specified shape 
and size | |4T]|122||143| . Eor neurofilaments, the OPD can be approximated by 
Guassian 1 122 143| . 

In previous experimental studies, the RDE and OPD of neurofilaments were calculated in 
selected regions of axonal cross-sections with almost no microtubules and organelles. To mimic 
such conditions, we performed simulations with exclusively neurofilaments, i.e., = nP = 0, 

and thus the only acting mechanisms are the pairwise repulsions and the Brownian motion of 
neurofilaments. We used a square domain with side length lp,m, and to minimize the effect of 
the boundary we used periodic boundary conditions. Under such conditions, the system Q 
reduces to 


dx 


N 


= T.K 

3,3 


NN . 
.3 


dt -f (Tn dwf 


1 < i < 


We initially put neurofilaments on a hexagon lattice inside the domain, and then ’’randomized” 
the distribution by simulating the model for sufficient time to observe no further change in the 
OPD or RDE To solve the model, we used the explicit Euler’s method with a time step 
h = 1/200 sec. 

We first investigated how the neurofilament distribution depends on its density. We took 
eNN = 0.5pN and used increasing neurofilament densities of 200 and 400 per fiuP (Eig. 
Rows A and B). Eor each case, the left panel is a plot of the coordinates of the neurofilaments 
after randomizing for 25 sec; the middle panel is a plot of the RDE which represent averages 
over 50 time frames between 25 sec and 30 sec; and the right panel is a plot of the averaged 
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Figure 3. Simulated neurofilament distributions with different densities and repulsion 
magnitudes. Left: snapshots of neurofilament positions after randomizing for 25 sec. Middle: 
the radial distribution functions (RDF, g{r)). Right: the bars are histograms of the occupancy 
probability distribution (OPD, pn) using randomly chosen circular windows with a radius 60 nm, 
and the black curves are their Gaussian fits. Middle and right plots represent averages over 50 
time frames between t = 25 sec and t = 30 sec. In all cases, e is short for . (A) = 200, 

e^N ^ Q 5 pN- (g) ^ 4QQ^ ^NN ^ Q 5 pjs^. ^ 4QQ^ ^NN ^ q_ 25 pN. other- 

parameters are the same as specified in Table [T] 


OPD and its Gaussian fit. The methods that we used to calculate the RDF and OPD are the same 
as in 1 122| and described in the supporting information (SI). These plots show that as the 
neurofilament density becomes higher, the separation of the peaks of the RDF becomes smaller, 
and the average and variance of the OPD becomes larger as the neurofilament density becomes 
larger. General features of these plots are in tight agreement with experimental data presented 
in ||4n [T22l[T43l . 


The magnitude of the repulsion between two neurofilaments depends on the charges of their 
sidearms. As the phosphorylation level of their sidearms becomes higher, their mutual repulsion 
becomes larger. We next investigated how the neurofilament distribution depends on the effect 
of sidearm phosphorylation by fixing the neurofilament density and varying . We took the 
neurofilament density to be 400 per pm^, and to be 0.25pN and 0.5pN. Fig. and C 
shows that as becomes larger, the locations of neurofilaments become more regular, the 
peaks of the RDF are better defined, and the variance of the OPD becomes smaller. 


Impairment of neurofilament transport leads to 
microtubule-neurofilament segregation. 


To investigate the mechanism of microtubule-neurofilament segregation in axons, we compared 
our simulations to experimental data obtained for IDPN in laboratory animals. We focused on 
IDPN because there is published data on both the rate and reversibility of the segregation. When 
IDPN is administered transiently by local injection into peripheral nerves, segregation appears 
within 2-6 hours and then disappears within 24 hours |^57 131| . 

Since neurofilament accumulation and axonal swelling occur on a much slower time course 


10 































A: normal axon 
t=1h 



B: 2h after IDPN 

t = 3h 



C: 5h after IDPN 
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Figure 4. Reversible segregation of microtubules and neurofilaments, neurofilament transport 
is blocked starting at f = 1 h and restored at f = 13 h. (A-F) Snapshots of the positions of 
microtubules, neurofilaments and organelles at f = 1 h, 3 h, 6 h, 12 h, 14 h, 20 h. All panels are 
from a single realization of the model. (A) microtubules and neurofilaments form a mixture 
under normal conditions. (B-D) blockage of neurofilament transport leads to gradual segregation 
of microtubules and neurofilaments. (E,F) restoration of neurofilament transport causes 
remixing of microtubules and neurofilaments. Small grey dots are free neurofilaments; small 
purple dots are neurofilaments engaged with microtubules; large black dots are microtubules; 
large blue circles are organelles. Parameters used: = 56, = 361. Neurofilament on-rate 

equals 0 between t = 1 h and 13 h. All other parameters are the same as in Table 


than the segregation, they can be ignored for the purposes of our current analysis. Therefore, for 
simplicity, we took to be a disk with fixed radius Rq — l/rm, and set the total number of 
neurofilaments to be constant. Specifically, if a neurofilament that was engaged with a 
microtubule left D, then it was replaced by a new neurofilament that entered D by association 
with a new randomly chosen microtubule. The total number of microtubules and neurofilaments 
in the domain were dete rmin ed based on the experimentally determined densities of 18//im^ and 
115/^m^, respectively jl3lj . We thus calculated by the formula = floor(187ri?Q) = 56 
and similarly we obtained = 356. Here the function floor(u) is the largest integer that is 
smaller than u. We considered organelles with b = 140 nm and a/b = 10 (Fig. |^). 

We started the simulations by including axonal transport of both neurofilaments and 
organelles, mimicking the conditions of normal axons. To distribute the neurofilaments and 
microtubules randomly without overlap, we first placed them on a hexagon lattice in D with no 
organelles, and then introduced volume exclusion and Brownian motion for enough time to 
randomize their positions. Starting from this initial condition, we then turned on the movement 
of both neurofilaments and organelles. Fig. is a snapshot of the simulated distribution of 
microtubules, neurofilaments, and organelles in a normal axon. The small grey dots are 
neurofilaments that are not engaged with microtubules, the small purple dots are neurofilaments 
that are engaged with microtubules, the large black dots are microtubules, and the large cyan 
circle is an organelle pushing into the cross-sectional domain. Note that a small fraction of the 
neurofilaments are bound to microtubules and moving along microtubules, that one microtubule 
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can transport multiple cargoes (neurofilaments or organelles), and that one organelle can engage 
with multiple microtubules simultaneously. 

We then blocked neurofilament transport selectively by resetting the binding rate of 
neurofilaments to microtubules, to be 0 at < = lb. This disengaged neurofilaments from 
their microtubule tracks and thus blocked their movement so that none could enter or leave D. 
Meanwhile, the transport of organelles was not affected; they continued to grab microtubules 
stochastically, pulling them together. This “zippering” effect caused the microtubules to 
gradually cluster (|^). By 6 hours, almost all the microtubules had migrated to the center of D 
and formed a single island surrounded by neurofilaments D). The central microtubule 
cluster contained organelles but relatively few neurofilaments whereas the peripheral zone of 
neurofilaments contained relatively few microtubules or organelles. This segregation pattern is 
strikingly similar to that observed in experiment and in disease, and the rate of segregation is 
comparable to that observed experimentally for local injection of IDPN into peripheral nerves of 
laboratory animals p4||131| . 

After observing segregation, we restored neurofilament transport by resetting to its 
original value at t = 13 h. This immediately allowed neurofilaments on the periphery of the 
microtubule core within a distance i?b of a microtubule to bind to that microtubule stochastically 
and then either unbind or exit D after a short while, as dictated by their stop-and-go transport 
behavior. As explained in the Methods, each neurofilament that exited D was replaced with a 
new neurofilament seeded adjacent to a randomly selected microtubule, but only if that 
microtubule was within a distance of Ri, from another neurofilament already in that plane. Over 
time this resulted in a gradual infiltration of neurofilaments into the microtubule cluster in a 
centripetal manner (i.e. from the outside edges progressing inward), leading to a gradual 
dispersal of the microtubules (|^) and a return to their normal interspersed organization Q^). 
These results agree tightly with previous experimental findings. 

To characterize the reversible segregation of microtubules and neurofilaments, we plotted the 
distribution and the mean of the pairwise distance between two microtubules (PDMT) as a 
function of time. Figs. and|^ are plots calculated from the simulation shown in Fig. 
which demonstrate a significant progressive decrease of the PDMT upon elimination of 
neurofilament transport {t = Ih) and subsequent increase upon restoration of neurofilament 
transport {t = 13h). Figs. HP an d|5p are plots for a normal axon for comparison. We see that 
under normal conditions, because microtubules and neurofilaments are interspersed, the 
distribution of PDMT is broad and the mean of it is about 0.8i?o; and as microtubules and 
neurofilaments segregate from each other, the distribution becomes more compact and the mean 
of the PDMT decreases by almost 40%. 

Another way to incorporate blockage of neurofilament transport is to increase the off-rate of 
neurofilaments We performed simulations with 100 times larger, and obtained 
similar results as in Fig. 

Microtubule zippering by moving organelles is the causal mechanism for 
segregation. 

In the above section we have shown that in the absence of neurofilament transport, organelle 
transport leads to microtubule-neurofilament segregation. As we noted earlier, organelles can 
interact with multiple microtubules simultaneously and thus pull or zip nearby microtubules 
closer together. We next investigated the importance of this zippering mechanism for the 
segregation of microtubules and neurofilaments. To do this, we introduced a maximum number 
of microtubules that a single organelle can interact with simultaneously, denoted by mmax, and 
investigated how the PDMT depends on Wmax in the absence of neurofilament transport. 

Fig. I^plots the mean of PDMT as a function of time given different values of rrimax- Each 
curve is averaged over five realizations with unpredictable seeds, and the error bars indicate the 
standard deviations over the realizations. If each organelle is only allowed to bind to one or two 
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Figure 5. Statistics of the pairwise distances between microtubules (PDMT). (A, B) IDPN 
treatment started at t=l hour and stopped at t=13 hours. (C, D) control. (A, C) distribution of the 
PDMT; data plotted for every 20 min. The pseudo color key represents the number of 
microtubule pairs. (B, D) mean of the PDMT; data plotted for every min. Parameters used are 
the same as in Fig. 


microtubules, i.e., Wmax = 1 or 2, then microtubules and neurofilaments remain mixed over 
time and segregation does not occur at all (blue and green). Indeed, for TOmax = 1, the mean of 
PDMT is slightly larger than that for a normal axon shown in Fig. |^. This is because organelles 
stir microtubules and neurofilaments and separate microtubules apart. As TOmax increases, the 
PDMT curve decreases faster and the time needed to reach complete segregation decreases. 
Scatter plots of microtubules and neurofilaments (not shown here) show that for TOmax = 4, 
partial but significant segregation was observed by 18 hours in all five realizations; for 
TOmax = 6, complete segregation was observed by 18 hours in four out of five realizations; for 
TOmax = 8 or 16, Complete segregation was observed in all realizations within 10 hours. These 
results suggest that microtubule zippering by moving organelles is the causal mechanism for the 
segregation of microtubules and neurofilaments in the absence of neurofilament transport. 

Dependence on the size and the flux rate of organelles. 

We next investigated the dependence of the segregation on the size and the flux rate of the 
organelles by simulating the model with different sizes of organelles, b = 140 nm or 6 = 70 nm, 
and different flux rates fc®. Fig. |^plots the mean PDMT over time for four situations: 
b = 140 nm and = 0.105/s (shown in blue), which is the same as in Fig. 5 = 140 nm 
and = 0.1575/s (shown in green); b = 70nm and = 0.105/s (shown in red); and 
b = 70 nm and = 0.21/s (shown in cyan). These results suggest that (1) for organelles of 
the same size, the more frequently they move through D, the faster the segregation occurs; (2) 
given the same flux rate across D, larger organelles are more capable of clustering microtubules 
and segregating them from neurofilaments than small organelles, and this is because on average 
larger organelles can interact with more microtubules simultaneously. 

Interestingly, simulations of the model demonstrate that during the segregation process 
microtubules frequently form smaller clusters first, then these small clusters gradually merge 
with each other to finally form a single large cluster near the center of the domain. These 
intermediate states were more apparent in simulations with small organelles, presumably 
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Figure 6. The effect of microtubule zippering by moving organelles. The mean of PDMT is 
plotted over time. The maximum number of microtubules that a single organelle can interact 
with simultaneously {rrimax) is set to be 1, 2, 4, 8, and 16 for the blue, green, red, cyan, purple, 
and yellow curves respectively. Each curve represents the average over 5 realizations of the 
model and the error bars are the standard deviation. All other parameters are the same as in Fig. 
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Figure 7. The segregation process depends on organelle size and flux rate. The mean PDMT is 
plotted over time. Organelle max radius b: 140 nm for blue and green curves (same as Fig. |^; 
b = 70 nm for red and cyan curves. Organelle flux rate same as Fig. |^(xl; blue and red), 
1.5-fold greater (xl.5; green) and 2-fold greater (x2; cyan). Each curve represents an average of 
5 realizations of the model and the error bars are the standard deviation. All other parameters are 
the same as in Fig. 


because the rate at which the smaller clusters merge is slower under this condition. Figs. [^-C 
are snapshots of these intermediate states captured in a single realization (cori'esponds to the 
cyan curve in Fig. |7]). A similar pattern of isolated clusters of microtubules has also been 
reported by Zhu et al (73) (Fig.[^) (see Discussion). 

Partial blockage of neurofilament transport: dosage effect. 

We finally investigated the cross-sectional distribution of microtubules and neurofilaments when 
neurofilament transport is partially blocked. In the case of segregation induced by IDPN, this 
might be considered equivalent to varying the IDPN concentration. To do this we reduced 
by different extents at f = Oh. Fig. |^plots the mean of PDMT over time for equals 0.5, 0.2 
and 0 times of its original value. Each curve was obtained by averaging over 5 realizations with 
unpredictable seeds, and the error bars indicate the standard deviations about the mean. The data 
indicate that when A;^ is small enough, there is insufficient neurofilament transport to counteract 
the organelle-dependent microtubule clustering, and segregation is observed. However, as A;^ 
becomes larger, the rate of microtubule clustering becomes slower and the resulting clusters 
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Figure 8. Intermediate states of the segregation process. (A-C) snapshots of the segregation 
process in a single realization. (A) microtubules form three clusters by f = 4h. (B) these 
clusters remain separated for several hours until two of them merge around t = 8h. (C) finally 
all microtubules form a single big cluster near the center of the domain. The dimension of the 
organelles: b = 70nm, ajb = 10. The flux rate of the organelles: 0.21/s. All other parameters 
are the same as in Fig. (D) multiple microtubule clusters observed in experiments with IDPN 
Big dots: microtubules; circular objects: organelles. Black arrows point to microtubule islands. 
Adapted from Fig. 9A of |[73|. The scale bar is 0.4 pm 


become less compact, reflecting less efficient segregation. Increasing has a similar effect: 
as k^j:j becomes larger, the rate of microtubule clustering becomes faster and the clusters 
become more compact (not shown). Thus the rate and extent of microtubule-neurofilament 
segregation is dependent on the extent of inhibition of neurofilament transport. 


Discussion 

Summary of our model. We developed a novel stochastic multiscale model for the 
cross-sectional distribution of microtubules and neurofilaments in axons. The model describes 
microtubules, neurofilaments, and membranous organelles as interacting particles in an axonal 
cross-section. It incorporates detailed descriptions of key molecular processes hat occur within 
seconds, including the axonal transport of neurofilaments and membranous organelles through 
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Figure 9. Dependence of microtubule-neurofilament segregation on Each curve plots the 
mean of the PDMT over time averaged over five realizations, and the error bars are the standard 
deviations. The rate is reduced to 50% (blue), 20% (green), and 0% (red) of the value in Fig. 
I^at t = Ih. All other parameters are the same as in Fig. 


this plane, as well as volume exclusion and Brownian motion of all the particles, and addresses 
the segregation phenomena that occur on a time scale of hours to days. The positions of the 
particles in the plane are governed by a system of stochastic differential equations. 

Mathematical models of the axonal transport of neurofilaments and organelles have been 
developed previously to describe the longitudinal distribution of cargoes along 
axons | fTT4HTT8l[T^[T^ . However, those models were in ID and did not consider the spatial 
arrangement and mechanical interactions of the cargoes and tracks in the radial dimension which 
are essential in understanding the segregation of microtubules and neurofilaments as well as the 
subsequent axonal swelling in neurological diseases. In our model, we describe in detail the 
dynamic interactions of neurofilaments, organelles, and nearby microtubules through molecular 
motors and volume exclusion in cross-section. Simulations of the model are in tight agreement 
with experimental data and generated a number of predictions that can be tested experimentally. 

Neurofilament and membranous organelle transports are competing processes. Simulations 
of the model demonstrate that if we block neurofilament transport selectively by preventing 
neurofilament binding to microtubules, while allowing organelle movement to continue, then the 
moving organelles tend to zipper nearby microtubules together so that they gradually segregate 
from the neurofilaments. The microtubule zippering action of the membranous organelles arises 
because we allow multiple motors to engage with a single organelle, which is consistent with 
experimental data and theoretical considerations |^90 121 129) . Restoration of neurofilament 
transport in the model allows the neurofilaments and microtubules to remix until their spatial 
distribution returns to normal. This suggests that neurofilament transport and organelle transport 
are competing processes in determining the cross-sectional distribution of microtubules: 
neurofilament transport can insert neurofilaments between adjacent microtubules, pushing those 
microtubules apart, while organelle transport can pull microtubules together when they move 
along multiple microtubules simultaneously, similar to a zipper. In normal axons, a dynamic 
balance between these two processes leads to the interspersed distribution of microtubules and 
neurofilaments, while in the absence of neurofilament transport, the microtubule zippering effect 
of organelle transport causes microtubules and neurofilaments to segregate. Thus our model 
predicts that the microtubule-neurofilament segregation that is observed in axons in neurotoxic 
and neurodegenerative diseases is a simple emergent property of the motile properties of 
membranous organelles that is triggered by selective impairment of neurofilament transport. An 
important and experimentally testable prediction of this study is that segregation is dependent on 
organelle movement. 

Why are the microtubule clusters mostly central? An intriguing feature of 
microtubule-neurofilament segregation, which is consistent across all published reports, is that 
the microtubules generally cluster in the center of the axon, surrounded by a peripheral band of 
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neurofilaments (see Introduction and Fig. [^. It is interesting to note that this was usually the 
case in our simulations also. According to our model, the segregation generated by microtubule 
clustering is caused by an exclusion of neurofilaments from the microtubule domain due to their 
failure to interact. The central location of the microtubule bundle is essentially a boundary effect 
which arises because microtubules at the periphery of the axon can only be pulled towards 
microtubules that are located more centrally whereas microtubules in the center can be pulled 
towards microtubules on all sides. The net result is that microtubule zippering by moving 
organelles tends to pull these polymers towards the axon center, displacing the neurofilaments to 
the periphery. The organelles co-segregate with the microtubules because they must follow the 
available tracks. 

Segregation proceeds via the merging of small microtubule clusters. An interesting 
observation in our simulations is that microtubule-neurofilament segregation tends to proceed 
initially via the formation of small microtubule clusters that subsequently merge together. This 
was more apparent in simulations with smaller organelles, which are less efficient at zippering 
microtubules together (see discussion below). Multiple small microtubule clusters have been 
reported in some studies on microtubule-neurofilament segregation induced by IDPN 154 ^ 
(see Fig. [^), but there is no published time course of segregation so it remains to be proven that 
these clusters are indeed intermediate states. Interestingly, microtubule zippering in our 
simulations also gives rise to the formation of small microtubule clusters in healthy axons. 
However, with ongoing neurofilament transport these clusters are transient and rarely merge to 
form larger ones. This is consistent with reports that small clusters of microtubules, often 
adjacent to one or more membranous organelles, are commonly observed in electron 
micrographs of axons 

Factors influencing rate of segregation. Our analysis gives us some insights into the factors 
that influence the rate of microtubule-neurofilament segregation. First, given the same number 
density, larger organelles are more effective at causing segregation, because they can interact 
with more microtubules simultaneously and they can pull together microtubules that are farther 
apart. Second, segregation occurs faster if the flux rate of the organelles is larger. Third, 
segregation occurs faster if the degree of neurofilament transport impairment is larger. These 
predictions are experimentally testable. It is also clear that there must be some dependence on 
the density of motors on the organelle surface, as well as the neurofilament:microtubule ratio. 
We are currently performing an extensive investigation on how the segregation phenomena 
depend on combinations of the model parameters using model simplification, 
nondimensionalization and mathematical analysis. These efforts will provide further insight of 
the biological problem and will be published elsewhere in the future. 

The predicted rate of segregation is comparable to that in real axons. The best experimental 
data on the kinetics of microtubule-neurofilament segregation is for animals treated with the 
neurotoxin IDPN. However, the rate of segregation in those animals depends on the mode of 
administration. When applied systemically to rats using a single intraperitoneal injection, 
segregation was first noted after 4 days, and after 4 such injections at 3 day intervals, the 
resulting segregation persisted for 6-16 weeks |m). In contrast, when applied locally at high 
concentration by sub-perineurial injection into peripheral nerve, microtubule-neurofilament 
segregation was evident after 2 hours, with the microtubule clusters becoming increasingly 
compact over the next 4-10 hours 124 Nagele et al. analyzed the pairwise distance 
between microtubules (PDMT) and observed full compaction by 8 hours after injection. Sixteen 
hours later, segregation was no longer seen in most axons, indicating an almost complete 
reversal | [5^ . In our simulations, we observed segregation within 4-12 hours of a complete 
cessation of neurofilament transport, and remixing within 2-8 hours after a complete resumption 
of neurofilament transport. This rate of segregation is comparable to the kinetics observed 
experimentally for injections of IDPN into nerves, and suggests that this delivery method results 
in a transient but acute inhibition of neurofilament transport. We predict that the slower time 
course of segregation that is observed when IDPN is administered systemically is due to the 
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lower effective dose experienced by the axons in those studies. The rate of remixing was a bit 
shorter in our simulations than in the experimental reports, which may be because we assumed 
an instantaneous recovery of neurofilament transport rather than a gradual one, which is more 
likely. 

What is the mechanism of neurofilament transport impairment? It is important to note that 
the impairment of neurofilament transport that leads to microtubule-neurofilament segregation in 
toxic neuropathies and neurodegenerative diseases also leads eventually to focal neurofilament 
accumulations and axonal swellings (see Introduction). Since microtubules are the tracks along 
which neurofilaments move, and since microtubule-neurofilament segregation appears early and 
precedes neurofilament accumulation and axonal swelling, it has been hypothesized that the 
segregation reflects the uncoupling of neurofilaments from their transport machinery 124 [28 |. 
Our modeling supports this hypothesis, but the molecular mechanism is unclear. Many of the 
neurotoxic agents that cause microtubule-neurofilament segregation and impair neurofilament 
transport (e.g. hexanedione, IDPN, carbon disulfide) are reactive molecules that could, or are 
known to, modify neurofilaments chemically |281. It is thought that these compounds react with 
specific amino acid residues to form protein adducts which may then modify protein interactions, 
and that such chemical modifications target neurofilaments preferentially or that they somehow 
render these polymers more susceptible than other cargoes to transport impairments. This 
selectivity could arise, for example, due to the unique structure or unusual amino acid 
composition of neurofilament proteins. The mechanism of impairment could be by interfering 
with their interaction with molecular motors or with the interaction of these motors with the 
microtubule tracks. Future experimental studies will be required to resolve such questions. 

How do neurofilament accumulations arise? The mechanism by which neurofilament 
accumulations arise is also of great interest given that this occurs in so many neurodegenerative 
diseases. Since local accumulations can only form if more neurofilaments move into a segment 
of axon than move out, the appearance of local swellings along axons implies some longitudinal 
instabilities in the transport of these cargoes. Therefore we propose that neurofilament 
segregation is an early event in neurofilament transport impairment but that longitudinal 
instabilities or non-uniformities in the transport impairment must arise to give rise to local 
accumulations and axonal swellings. We plan to address this in future studies. Due the complex 
spatial and temporal nature of this problem, which entails the interactions of multiple dynamic 
components, we believe that a full understanding can only be achieved by a combination of 
experimental and modeling approaches. Our present study is an important first step. 
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Supporting Information A: calculation methods for the RDF 
and OPD in Fig. 3 

To calculate the RDF, for each time frame we used the centers of the neurofilaments in the 
domain [—0.2/rm, 0.2/rm] x [—0.2/im, 0.2p,m] as reference points. We chose to use reference 
points in a smaller domain in order to avoid boundary effects. For each reference point we 
binned the center-to-center distances (r) between the reference neurofilament and all other 
neurofilaments with a bin size Ar = Inm, and normalized by the factor 27rr where is 
the density of neurofilaments in the whole domain. We then took the average over all reference 
points for each time frame and then over all time frames to obtain g{r). We note that g{r) is 
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noisy due to under sampling and the small bin size (1 nm). Averaging over more time frames or 
using a larger bin size gave s moo ther g{r). To compute the OPD pn, we used circular windows 
with fixed radii 60 nm as in 11221. For each time frame, we sampled centers of lOn^ circles 
according to the uniform distribution in the domain [—O.S/rm, 0.3/rm] x [—0.3/im, 0.3/im] and 
calculated the particle occupancy number for each circle. We then produced a histogram of 
occupancy numbers obtained for all circles and normalized it by the total number of circles to 
obtain 


Supporting Information B: parameter estimation 

The parameters of the model are summarized in Table 1. In the following, we explain the 
methods used to estimate the parameters. 


Kinetic rates for neurofilament and organelle transports. 

We assume that a neurofilament can only bind to a single microtubule at one time whereas an 
organelle, which is much larger, can interact simultaneously with one or more microtubules. The 
capturing radius Rf, for these interactions is taken to be 80 nm, which is the length of a kinesin 
motor when fully extended, as measured by electron microscopy 11231 . 

The movement of neurofilaments along axons has been modeled as a bidirectional 
independent velocity jump process in ID with two distinct pausing states in which a 
neurofilament either moves anterogradely or retrogradely, pauses for a short time, or pauses for a 
long time | |124||125) . Since that model is one dimensional, the kinetic rates extracted represent 
space averages over all neurofilaments in the axonal cross-section. For this reason, the kinetic 
rates for neurofilament transport in our model are related, but not identical, to the rates extracted 
in those other studies. We estimate the rate at which a neurofilament binds to a nearby 
microtubule within distance Ri, in our model, to be five times as la rge a s the transition rate 701 , 
the rate for a short-pause neurofilament to start moving, defined in ||l25l. We estimate k^ff, the 
unbinding rate of neurofilaments in our model, to be th e sam e as the rate 710 , the rate for a 
moving neurofilament to stop and pause, as defined in 1 125) . Using 701 = 2.0 x 10“^/s and 
710 = 6.5 X 10“^/s for myelinated axons, we obtain = 1.0 x 10“^/s and 
k^jj: = 6.5 X 10“^/s. 

The rate that a neurofilament leaves the domain, is estimated in the following way. The 
average time for a neurofilament to move through D is , where is the average length 

of moving neurofilaments and is the speed of the filament. Assuming that the neurofilament 
departure events are exponentially distributed then the rate is the reciprocal of the average 
time, i.e., k^^^ = /Lm- Taking Lm « 5/im, which is the appr oximate a verage length of 
moving neurofilaments in cultured neurons, and « 0.5/im/s 1126 1271, we obtain 


-1 


=0-ls 

The binding and unbinding rates of organelles to a nearby microtubule, and k^^p are 
assumed to be constant 2 /s which are comparable to the rates used in previous mathematical 
models of vesicular transport along microtubules [128 129) . 

The passage rate for the organelles is calculated using their cross-sectional density. 
Assuming the cross-sectional organelle density is then the total number of organelles in D is 
given by ttRqP^ . Denoting the speed of the organelle along microtubule to be s'^ and the length 
of the organelle to be 2a, then the time that an organelle remains in D is given by 2al. 
Assuming that organelle arrival is a Poisson process, then the Poisson rate can be calculated as 
ttRqP^ divided by 2als^, that is, 

k?n = 7ri?oP'^s°/(2a). 

The sizes and densities of axonal organelles have been most carefully studied in ]130[ . 
According to Table 1 in 11301, there are two major sizes of organelles: large ones such as 
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mitochondria have an average cross-sectional diameter of 280 nm, and small ones classified as 
tubular and vesicular profiles have an average cross-sectional diameter of 50 nm. The densities 
of these organelles were counted in longitudinal slices of axons, and need to be converted to 
cross-sectional densities for our model. We note that the cross-sectional density and the 
longitudinal density (denote as pf) are generally different. However, the area fraction occupied 
by the organelles averaged over cross sections and longitudinal sections are comparable, and this 
relation can be used to convert the longitudinal density to cross-sectional density. Assuming that 
the objects are cylinders with radius r and length I, the area fractions measured in cross and 
longitudinal sections are given by p‘^ and 2rlpf respectively. Equating these two, we obtain 
p^ = 2lpf’ /(Trr). For the large organelles, 2r = 280nm, I = 870nm, pf = 0.023/pm^, and 
this leads to p^ = 0.091//im^. For the small organelles, 2r = 50nm, I = 180nm, 
pP = 0.02/pm^, and this leads to p^ = 0.0917/pm^. 

The density of organelles estimated above is comparable to the density data in published 
IDPN studies. For example, in Table 1 of | |13H , the mitochondrial density was measured to be 
0.187 — 0.250/pm^. Taking p^ = 0.187/p,m , a = 10& = 2.8p,m, = lp,m/s, and 

i?o = we obtain fc® = 7ri?Qp‘^s‘^/(2a) = 0.105/s. 


Parameters in the pairwise repulsion and elastic spring forces. 


We first estimate Lr which is the maximum distance for pairwise repulsions of particles. The 
parameter Lr for neurofilament-neurofilament repulsion is approximately 21^, where is the 


equilibrium brush thickness of the neurofilament sidearms. This was given by Eqn. (4) in 11221, 
that is, 

'12y/^_. /a\2/3 

.2 ) 
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Here N is the number of amino acids per sidearm, a is the length of an amino acid and s is the 
spacing betwe en neurofi lament sidearms. Taking N = 679 | |122) , a = 3.5 angstroms | |132||133| , 
and s = 3 nm |l22 1341, we obtain « 60.6 nm and thus Lr = 121.2 nm. For simplicity, we 
use the same Lr for pairwise interactions of all kinds of particles and their interactions with the 
domain boundary. 

We estimate the force prefactors in the following way. We assume that the repulsion force 
between two neurofilaments is approximately 1 pN when their surface distance is d = 40 nm. 
Under this assumption we have = 1/ {Lr/d — 1). Taking Lr = 121.2 nm, we obtain 
^NN ^ Q g denote by Sr for simplicity of notation. We assume that the force 

prefactor for microtubule-microtubule, microtubule-neurofilament repulsions and repulsions of 
microtubules and neurofilaments with the boundary are the same as e^. For repulsions that 
involve organelles, we use a prefactor that is five times as large, that is, = 5£r- 

Organelle movement can cause significant flow of the axoplasm near their surfaces and 
displace nearby microtubules and neurofilaments. As an organelle pushes into D, its radius 
increases and it pushes nearby fluid and particles away from itself; as it moves away from D, 
instead of leaving void behind it, it creates negative pressure which draws the axoplasm to flow 
back and fill the space. In this model, we do not include the hydrodynamic interactions among 
these particles explicitly, but include this effect by adjusting the force factors e associated with 
organelles. Specifically, when organelles push into the domain, we double and to take 
into account the contribution of the fluid flow. 

The effective spring constants and are calculated in the following way. The mean 
surface distance between a microtubule and a cargo engaged on it has been observed to be 17 nm 
(denote as Iq) \ 135 136) . We assume that the spring force and the repulsive force between a 
microtubule-neurofilament pair equilibrates at Iq, i.e., k^Iq = £r{Lr/lo — 1). From this 
assumption we obtain 




— ^r(,Lr/I q 1)/^0- 


(S2) 
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Plugging the values of £r, and /q into this expression, we obtain = 0.18pN/nm. 
Similarly, we assume that the spring force and the repulsive force of a microtubule-organelle 
pair equilibrates at Iq, and this leads to 

K° = 5e,(L,Ao-l)Ao, (S3) 


which is /t® =0.9 pN/nm. We note that the spring constant for a single kinesin motor is 
estimated to be 0.2 — 0.4 pN/nm in 11371. The spring constants here are different from the 
spring constant of a single motor in two ways: first, when a cargo moves along a microtubule 
there could be multiple motors being active, and thus the spring constants used here represent 
the sum of the spring constants over all active motors; second, the spring constants used here 
only take into account projections of the elasticity of individual motors in the plane orthogonal 
to the microtubule, whereas molecular motors are most likely slanted when dragging a cargo 
along a microtubule. 


Drag and diffusion coefficients 


We treat neurofilaments as slen der c ylinders, and estimate the drag coefficient per unit length 
by the formula given in 11381, = 47 rp/(ln(LAr/ 6 ) + In 2 — 0.5), where 77 is the 


viscosity of the axoplasm, is the characteristic length of the cylinder, and b is the 
characteristic radius of the cross-section. Th e viscosity of the cell cytoplasm is estimated to be 
p ^ 3 — 5 pN • s/pm^ for mammalian cells 1139 


140| . The persistence length of 


neurofilaments is 200-450 nm | |141||142) . The characteristic radius of a neurofilament with its 
sidearms is about 20 nm p43[ . Taking Tat = 500nm , b = 20nm, we obtain 

/L = 47rp/(ln(Ljv/6) -f ln2 — 0.5) « 14.7 pN • s/ Taking the length of the 
neurofilament to be L = 5/im, we obtain « 73.5 pN • s/ fim. 

We estimate the drag coefficient for microtubule in a similar way. The persistence length of 
microtubules has been estimated to be 80 ± 20 /im 044). We estimate the characteristic radius 
of the microtubule with its associated proteins in axons to be 37.5 nm. Taking 
rj = 4pN • Lm = 80/rm and b = 37.5 nm, we obtain = 6.4pN • s//rm^. We 

note that this estimate is close to the estimate obtained by treating microtubules as infinitely long 
cylinders, for which one can use the Oseen drag formula /L = 47177 / log(4/i?e — 7 + 0.5) 
where 7 « 0.5772 is the Euler’s constant and Re is the Reynolds number P45[|146| . The 
characteristic flow velocity in axonal cross-section is approximately U = 0.2/im/s given by the 
pushing of the organelles. Assuming that the density of the cytoplasm p to be the same as water, 
we find that Re = 2Ubp/r] « 4.3 x 10“^. Plugging rj and Re into the Oseen’s formula, we 
obtain the drag per unit length to be 6.3 pN • s/ pm^. Taking the sectional length of microtubule 
to be L = 80/rm, we obtain « 6.4 x 80 = 512 pN • s/ pui. 

We estimate the drag coefficient for organelles using the formulas for a prolate ellipsoid 
given in P46) . Assume that the major axis of the ellipsoid is 2a and the minor axis is 2b, then 
the drag coefficient per unit length is given by 


p^/L = 




for 0 < e < 1 , 


2e -I- (Se^ - 1) ln[(l -f e)/(l - e)] 

where e = y^(1 — (b/a)^) is the eccentricity. For b/a = 0.1, 0.2 and 0.5, we have 

p^/L « 14.4, 17.9 and 26.0 pN • s/pm^. For organelles with b = 140nm and b/a = 0.1, we 

have p^ « 40.3pN • sec/pm. 

We calculate the diffusion coefficients with k = M, N, and O using the Einstein relation 


Dk = ksT/pk, 


(S4) 


where fcs is the Boltzmann’s constant and T is the absolute temperature. At room temperature 
(25 °C or 298 K), one has = 4.11pN • nm. Using pjq = 73.5 pNs / pra, we obtain 
Dm « 5.6 X 10“® pm^/s. We calculated Dm and Dq in a similar way. Finally cr^ is given by 
the relation ak = \/2JJk. 
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Supporting Information C: simulation algorithm 


To solve the model numerically, we chose a time step h much smaller than all the time scales 
involved in Mechanisms 1-3, treated the binding and unbinding, arrival and departure of cargoes 
explicitly at discrete time steps, and integrated the model system (Eqn. 8 in main text) using the 
explicit Euler’s method. Because cr^, k = M, N, C are constant in time, the numerical 
integrator has strong order 1.0 1 147) . Eor the segregation simulations, we used h = 1 /50 sec if 
there was no organelle in D\ otherwise we used h = 1/1600 sec in order to deal with the 
stiffness of the equations introduced by the pushing of organelles when they move into D. The 
algorithm for a typical time step is summarized below. 


Algorithm for the model 

1. Stochastic removal of moving neurofilaments. Eor each moving neurofilament, generate a 

random number r uniformly distributed in [0, 1]. If r < 1 — then remove it from 

D. 

2. Update zf and r® for each organelle in D according to Eqn. 1 and Eqn. 4 in the main 
text. If zf becomes bigger than or equal to a then remove the /-th organelle and release all 
microtubules from it. 

3. Stochastic unbinding of cargoes, i.e., neurofilaments and organelles, from their engaged 

microtubules. If a microtubule and an organelle are engaged, then generate a random 
number r uniformly distributed in [0 1], and let them unbind if r < 1 — They 

also unbind if their surface distance is bigger than Ri,. The same method is used to update 
microtubule-neurofilament binding. 

4. Stochastic binding of cargoes to microtubules. 

(a) If the surface distance of a cargo and a microtubule is smaller than i?^, then generate 
a random number and determine if they intend to bind to each other. Loop through 
all microtubule-neurofilament and microtubule-organelle pairs, and find all potential 
binding events. 

(b) Accept or reject the potential binding events according to the availability of the 
associated microtubules in a random order. We assume that one microtubule has 5 
tracks, each neurofilament occupies one track, an organelle with maximum radius 
140nm occupies 2 tracks, and an organelle with maximum radius 70nm occupies 
1.5 tracks. 

5. Addition of new neurofilaments and organelles to D. The number of new neurofilaments 

equals the number that has been removed in step 1. The number of new organelles is 
determined using the rate Since the time step h is much smaller than the mean arrival 
time of organelles we introduce a new organelle in each time step with probability 
(1 — We then add these cargoes at Iq = 17 nm away from a randomly chosen 

microtubule at a random angle. If the random location overlaps with an existing particle or 
the associated microtubule is too crowded, then a different microtubule and an angle is 
generated randomly. Once a new cargo is added, it is bound to the selected microtubule. 

6. Update the positions of the microtubules, organelles and neurofilaments in D by 
integrating the model system (Eqn. 8 in the main text) using the explicit Euler’s method, 
i.e., 

x/(f + h) = x/(<) -f -f CTfcrfv^, 

l<i<n^, k = M,N,F. 
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where are pairs of random nu mber s generated from the standard normal distribution 
using the Box-Muller transform jl48) . To avoid large values of y\, they are regenerated if 
the absolute value of any component is greater than 5. 

7. Go to Step 1 for the next time step. 
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